//
//
//

#include "simpleMatrixMulter.h"

#include <stdio.h>

#include <cutil_inline.h>

extern bool verbose;

// Matrix multiplication - Host code 
// Matrix dimensions are assumed to be multiples of BLOCK_SIZE 
void MatMul(const Matrix A, const Matrix B, Matrix C) 
{ 
	// Load A and B to device memory 
	Matrix d_A; 
	d_A.width = A.width; 
	d_A.height = A.height; 
	size_t size = A.width * A.height * sizeof(float); 
	cudaMalloc(&d_A.elements, size); 
	cudaMemcpy(d_A.elements, A.elements, size, cudaMemcpyHostToDevice); 
	
	Matrix d_B; 
	d_B.width = B.width; 
	d_B.height = B.height; 
	size = B.width * B.height * sizeof(float); 
	cudaMalloc(&d_B.elements, size); 
	cudaMemcpy(d_B.elements, B.elements, size, cudaMemcpyHostToDevice);

	// Allocate C in device memory 
	Matrix d_C; 
	d_C.width = C.width; 
	d_C.height = C.height; 
	size = C.width * C.height * sizeof(float); 
	cudaMalloc(&d_C.elements, size); 
	// Invoke kernel 
	dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE); 
	dim3 dimGrid(B.width / dimBlock.x, A.height / dimBlock.y); 
	
	MatMulKernel<<<dimGrid, dimBlock>>>(d_A, d_B, d_C); 
	// Read C from device memory 
	cudaMemcpy(C.elements, d_C.elements, size, cudaMemcpyDeviceToHost); 
	// Free device memory 
	cudaFree(d_A.elements);
	cudaFree(d_B.elements);
	cudaFree(d_C.elements);
} 

// Matrix multiplication kernel called by MatMul() 
__global__ void MatMulKernel(Matrix A, Matrix B, Matrix C) 
{ 
	// Each thread computes one element of C 
	// by accumulating results into Cvalue 
	
	float Cvalue = 0; 
	int row = blockIdx.y * blockDim.y + threadIdx.y; 
	int col = blockIdx.x * blockDim.x + threadIdx.x; 
	for (int e = 0; e < A.width; ++e) 
		Cvalue += A.elements[row * A.width + e] * B.elements[e * B.width + col]; 
	C.elements[row * C.width + col] = Cvalue; 
}

static void initMatrixSequence (Matrix &A)
{
	A.width = A.height = BLOCK_SIZE * 64;
	A.elements = static_cast<float *>(malloc(A.width * A.height * sizeof(*A.elements)));

	for(int i = 0 ; i < A.width * A.height ; ++i)
		A.elements[i] = i + 1;
}


void dumpMatrix(Matrix &A)
{
	printf("Matrix (%d x %d)\n", A.height, A.width);
	for(int i = 0 ; i < A.height ; ++i) {
		for (int j = 0 ; j < A.width ; ++j) {
			printf("%6g ", A(i, j));
		}
		printf("\n");
	}
}

void testMatMul()
{
	Matrix A;
	initMatrixSequence(A);
	if(verbose)
		dumpMatrix(A);

	Matrix B;
	initMatrixSequence(B);
	if(verbose)
		dumpMatrix(B);

	Matrix C; 
	C.width = A.width; 
	C.height = A.height;
	C.elements = static_cast<float *>(malloc(C.width * C.height * sizeof(*C.elements)));

	unsigned int timer = 0;
	cutilCheckError(cutCreateTimer(&timer));
	cutilCheckError(cutStartTimer(timer));
	MatMul(A, B, C); // ************* RUN THE TEST *************
	cutilCheckError(cutStopTimer(timer));

	if(verbose)
		dumpMatrix(C);

	printf("C[%d, %d] == %g\n", C.height - 1, C.width - 1, C(C.height - 1, C.width - 1));
	printf("Processing time: %f (ms)\n", cutGetTimerValue( timer));
	cutilCheckError(cutDeleteTimer(timer));
}
